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Abstract 

We  describe  a  method  based  on  the  Total  deviation  approach  whereby  we  improve  the  confidence  of 
the  estimation  of  the  Hadamard  deviation  that  is  used  primarily  in  GPS  operations .  The  Hadamard- 
total  deviation  described  in  this  paper  provides  a  significant  improvement  in  confidence  indicated  by 
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an  increase  of  1.3  to  3.4  times  the  one  degree  of  freedom  of  the  plain  Hadamard  deviation  at  the 
longest  averaging  time.  The  new  Hadamard-total  deviation  is  slightly  negatively  biased  with  respect 
to  the  usual  Hadamard  deviation ,  and  r  values  are  restricted  to  less  than  or  equal  to  T/ 3,  to  be 
consistent  with  the  usual  Hadamard’s  definition.  We  give  a  method  of  automatically  removing  bias 
by  a  power-law  detection  scheme.  We  review  the  relationship  between  Kalman  filter  parameters  and 
the  Hadamard  and  Allan  variances ,  illustrate  the  operational  problems  associated  with  estimating 
these  parametersy  and  discuss  how  the  Hadamard-total  variance  can  improve  management  of  present 
and  future  GPS  satellite  clocks. 


1  INTRODUCTION 

Using  a  type  of  Hadamard  variance,  the  goal  of  this  paper  is  to  reduce  the  uncertainty  of  long-term 
estimates  of  frequency  stability  without  increasing  the  length  of  a  data  run.  For  measurements  of 
frequency  stability,  the  two-sample  frequency  variance  known  as  the  Allan  variance  was  generalized 
to  an  iV-sample  variance  weighted  with  binomial  coefficients  by  R.  A.  Baugh  [1].  The  case  of  the 
three-sample  frequency  variance  that  is  used  here  is  the  Picinbono  variance  [2]  times  |.  However, 
in  this  paper,  it  will  be  called  a  Hadamard  variance  (following  Baugh’s  work)  that  is  defined  as 
follows.  Given  a  finite  sequence  of  frequency  deviates  {yn,n  =  1, . . . ,  Nyman,},  presumed  to  be  the 
measured  part  of  a  longer  noise  sequence  and  with  a  sampling  period  between  adjacent  observations 
given  by  To,  define  the  r  =  mro-average  frequency  deviate  as 

^  m—1 

Pn(m)  =  —  Vn+r  w 

m  3=0 

Let  Hn(m)  —  yn (rn)  —  2yn+rn(m) + Vn+2m (m)  be  the  second  difference  of  the  time-averaged  frequen¬ 
cies  over  three  successive  and  adjacent  time  intervals  of  length  r.  Define  the  Hadamard  variance 
as 

H%2(r)  =  (2) 

where  <  •  >  denotes  an  infinite  time  average  over  n,  and  Hffy2  depends  on  m. 

The  GPS  program  office  uses  this  particular  time-series  statistic  for  estimating  Kalman  al¬ 
gorithm  coefficients  according  to  [3],  which  coefficients  will  be  discussed  in  a  later  section.  The 
Hadamard  deviation  H<Ty(r)  a  function  that  can  be  interpreted  like  the  more  efficient  Allan  de¬ 
viation  as  a  frequency  instability  vs.  averaging  time  r  for  a  range  of  frequency  noises  that  cause 
different  slopes  on  }\ay  (r).  This  is  shown  in  figure  1.  For  estimating  Kalman  drift  noise  coefficients, 
h o-y(r)  is  inherently  insensitive  to  linear  frequency  drift  and  reports  a  residual  “noise  on  drift”  as 

3 

a  T2  slope,  or  what  is  commonly  called  random  run  FM  (RRFM).  This  is  in  contrast  to  the  Allan 
deviation,  which  is  sensitive  to  drift  and  causes  a  r+1  slope.  If  the  level  of  drift  is  relatively  high, 
it  masks  the  underlying  random  noise.  It  is  customary  to  estimate  and  remove  overall  frequency 
drift.  Depending  on  the  method  of  drift  removal,  this  procedure  can  significantly  alter  the  Allan 
deviation  in  the  longest  term  r  region  of  interest,  so  estimating  underlying  noise  can  be  a  formi¬ 
dable  task  for  any  given  data  span.  On  the  other  hand,  the  Hadamard  deviation  is  unaffected  by 
removing  overall  frequency  drift.  For  this  reason,  it  is  the  preferred  statistic  in  situations  in  which 
the  frequency  drift  may  be  above  the  random  noise  effects,  which  is  the  case  with  the  use  of  Rb 
clocks  in  the  GPS  Block  II  satellite  program.  We  do  not  imply  that  systematics  such  as  frequency 
drift  can  be  ignored.  Indeed,  satellite  clocks  are  changed  and  these  systematics  must  be  learned  as 
quickly  as  possible  to  ensure  a  smooth  changeover. 
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FREQUENCY  STABILITY 


Sample  Time,  %,  Seconds 


Figure  1:  The  Hadamarcl  deviation  (root  Hvar)  shows  FM  power-law  noises  as  straight  lines  in 
addition  to  PM  sources  of  noise  for  r-domain  power-law  exponent  /x  (that  is,  h ay(T)  °c  tm)  range 
of  —  2  <  /x  <  3.  We  define  a  new  estimator  that  can  be  interpreted  identically  called  Hadamard- 
total  deviation  (root  TotHvar)  and  that  has  significantly  improved  confidence  at  long  term.  The 
Hadamard-total  deviation  is  insensitive  to  linear  frequency  drift  that  can  mask  characteristic  ran¬ 
dom  noise  typically  encountered  here  in  the  region  where  r  =  one-week  and  longer.  The  goal  is  to 
identify  /x  even-integer  power-law  noises  and  accurately  estimate  their  levels  in  order  to  set  system 
parameters  associated  with  the  GPS  Kalman  filter. 


Throughout  this  writing,  we  will  make  comparisons  using  the  traditional  best  statistical  estima¬ 
tors,  denoted  by  “Hvar”  and  “Avar”  referring  to  the  maximum-overlap  estimators  of  the  Hadamard 
and  Allan  variances.  Section  2  reviews  the  “total”  approach  to  improving  statistical  estimation. 
Sections  3  and  4  give  two  methods  of  computing  total  Hadamard  variance,  designated  as  TotHvar, 
using  measurements  first  of  fractional  frequency  deviations  and  then  of  time  deviations.  Then  we 
quantify  the  advantage  of  TotHvar  over  Hvar  in  Section  5,  giving  formulae  for  computing  bias  and 
equivalent  degrees  of  freedom  (edf)  of  TotHvar.  Section  6  gives  a  method  for  efficiently  determining 
the  noise  type  at  a  given  r-value  for  automatically  correcting  the  bias  and  determining  confidence 
intervals  for  the  range  of  noises  considered  by  TotHvar.  Section  7  reviews  how  an  estimate  of 
r-domain  frequency  stability  is  used  to  set  Kalman  filter  parameters  (or  q’ s)  used  in  GPS  opera¬ 
tions,  problems  associated  with  the  application  of  either  the  traditional  Allan  variance  or  Hvar  to 
the  Kalman,  filter,  and  how  TotHvar  serves  as  a  unifying  solution.  Finally,  Section  8  discusses  a 
past  scenario  in  GPS  operations  in  which  TotHvar  is  applied  to  real  data  showing  the  benefit  of 
improved  estimation  of  long-term  frequency  stability. 

2  THE  “TOTAL”  APPROACH 

The  total  estimator  approach  has  been  developed  to  improve  confidence  of  major  statistical  tools 
used  in  analyzing  and  characterizing  instabilities  in  phase  and  frequency  of  oscillators  and  syn¬ 
chronization  systems  [4-9].  Making  a  “total”  estimator  of  eqn.  (2)  involves  joining  each  real  data 
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subsequence,  namely  the  subsequence  of  iji  that  goes  into  each  Hn(m )  term,  at  both  its  endpoints 
by  the  same  original  data  subsequence  so  that  it  repeats.  This  creates  a  new  extended  version  of 
each  yi  subsequence  that  may  be  extended  by  a  forward  or  backward  repetition,  with  or  without 
sign  inversion,  thus  with  four  possible  ways  to  extend.  From  numerous  simulation  studies,  we  have 
determined  that  an  extension  by  even  (uninverted)  mirror  reflection  of  linear-frequency-detrended 
Hn{m)  subsequences  yields  the  largest  edf  gain  and  least  bias  for  the  range  of  noise  types  identified 
by  standard  Hvar. 


3  COMPUTATION  USING  ?/, ,-SERIES 

Hn(m)  is  computed  from  a  3m-point  data  segment  or  subsequence  {'Ui}n  =  {yi, i  =  n, . . .  ,n  +  3 m  — 
1}.  Before  applying  any  data  extensions,  we  must  remove  a  linear  frequency  trend  (drift)  from  each 
subsequence  by  making 

°Vi  =  Vi-  c\i , 

where  c\  is  a  frequency  offset  that  is  removed  to  minimize  —  °Z/*)2>  to  satisfy  a  least- 

squared-error  criterion  for  the  subsequence.  In  practice,  it  is  sufficient  to  compute  this  background 
linear  frequency  slope  by  averaging  the  first  and  last  halves  of  the  subsequence  divided  by  half 
the  interval  and  subsequently  subtracting  the  value.  Now  extend  the  “drift-removed”  subsequence 
at  both  ends  by  an  uninverted,  even  reflection, 
extensions  as  follows.  For  1  <  l  <  3m,  let 

°yt  i  =  °yn+i~ i,  °yi+3m+i- 

to  form  a  new  data  subsequence  denoted  as  {°yf}n  consisting  of  the  drift-removed  data  in  its  center 
portion,  plus  the  two  extensions,  and  thus  having  a  tripled  range  of  n  — 3m  <i<  n+6m— 1  with  9m 
points.  To  be  clear,  we  now  have  extended  subsequence  {°yf  }n  ~  {°yf ,  i  =  n— 3m, . . . ,  n+6m— 1}. 


Utility  index  l  serves  to  construct  the 

Vn+Qm — l >  v*0 


Define 


Total]  i<7y(m,  T(),  btymax) 


^  Nymax  1 

_ L_ _  v 

6  (Nymax-  3m+l) 


1  n-\-3m—  1  9\ 


6m  . 

v,  z—n—6m 


(4) 


for  1  <  m  < 


Nv 


,  where  [c\  means  the  integer  part  of  c  and  notation  °Hf  (m)  means  that 

Hn(m)  above  is  derived  from  the  new  triply-extended  subsequence  \°yf}-  The  symmetries  of  the 
extension  and  the  Hvar  filter  allow  the  computational  effort  to  be  halved.  Let  k  =  [3m/2j.  We 
need  to  calculate  only  for  n~k<i<n+k+  3 m  —  1,  and  (m)  only  for  n  —  k<i  <~n  +  k. 

Then 


n+3m-l 


n-\~k—  1 


— i  2  to  \  ru  4.  2  /  \  2  /  \  ^ 

J2  Hi  (m))  =  2  £  (#H?(m))  +  (m)J  +  [*H°+k(m)J  ,m  even, 

i—n—k-\r  1 

ii+fc  2 

=  ^  £  (**?("»>)  ,m  odd.  (5) 


i—n—k 
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4  COMPUTATION  USING  x„-SERIES 


The  methodology  described  above  can  be  written  in  terms  of  calculations  on  residual  time  differences 
between  clocks,  namely  an  Xi-series  (to  adhere  to  usual  notation),  recalling  that 

Vi  (m)  =  (xi+m  -  Xi)  /  (mr0) . 

Thus  in  the  total  approach  applied  to  x^-series,  the  data  extensions  on  subsequences  of  x*  will  be 
constructed  in  such  a  way  that 

°vt  =  ("4.  -°4)  /TO, 

in  agreement  with  section  3  above.  This  has  the  effect  of  requiring  an  odd  mirror  extension  and  a 
third-difference  operator  when  considering  subsequences  of  Xj.  The  Hadamard  variance  discussed  in 
section  3  as  a  second-difference  operator  on  r-averaged  yn  values  can  now  be  re-expressed  in  terms 
of  a  third-difference  operator  on  time-error  x*- values.  The  sample  variance  (or  mean  square)  of 
these  third  differences  falls  neatly  into  a  class  of  structure  functions,  namely  the  variance  produced 
by  a  difference  operator  of  order  three  [10].  The  modified  Allan  variance  can  also  be  treated  as  a 
third-difference  variance  [11]. 

The  Xj-subsequence  that  corresponds  to  the  y^-subsequence  starting  at  n  is  {xj,  n  <  i  <  n  +  3m} , 
which  has  3m  +  1  terms.  Compute  the  detrended  subsequence  °Xi  according  to 

Xn  /‘ri  \  k  X,H  3m— k  T  X,i-|-3m. 

k  (3m  —  k )  ’ 

°Xi  =  Xi  —  \c2  [i  —  n)(i  —  n  —  3m) ,  n<i<n  +  3m. 

Define  the  extended  subsequence  j  °xf,  n  —  3m  <i<n  +  6m  j  by 

°xf  =  °Xi,  n  <  i  <  n  +  3m, 

°xt_l  =  2  (°xn)  -  °xn+i ,  1  <  l  <  3m, 

^'It.+Snn+l  =  2  (  X„  \-Srii)  ®n+3m- h  1  —  ^  —  3m. 

Then 

mr0  (° H*  (m))  =  -  °xf  +  3  (°sf+m)  -  3  (°x#  2m)  +  °xf+3m,  n  -  3m  <  i  <  n  +  3m  -  1, 

where  °Hf  (m)  has  the  same  meaning  as  in  Section  3.  Now  the  Hadamard-total  variance  is  com¬ 
puted  from  (4)  as  before  with  Nymaas  =  Nxm(lx  —  1.  Because  of  symmetry  we  need  #x?  only  for 
n  —  k<i<n  +  k  +  3m,  and  (5)  applies. 

5  BIAS  AND  EQUIVALENT  DEGREES  OF  FREEDOM 

We  consider  the  random  frequency-modulation  (FM)  noises  since  these  dominate  at  long-term 
averaging  times  where  we  can  capitalize  on  the  improved  confidence  of  using  the  total  approach. 
To  analyze  phase-modulation  (PM)  noises,  one  would  usually  use  Total  TDBV  [6]  rather  than 
the  Hadamard  deviation.  For  brevity,  let  Total  a  ^(m,  tq,  Nyrnax)  be  TotHvar  (r,  T),  where  r  = 
mro,  T  =  Nyrnaxro.  The  normalized  bias  and  edf  for  TotHvar  are  given  by 
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nbias (r)  = 


£  {TotHvar  (r,T)} 
L  E  {Hvar  (r,  T)} 


-  —  1  =  a, 

(6) 

T/t 

bo  +  bir/T 

(7) 

where  E{-}  is  expectation  of  {•},  0  <  r  <  r  >  16ro  (to  be  explained),  and  a,  bo,  and  b\  are  given 
in  Table  1  for  the  five  FM  noise  types  considered  by  the  Hadamard  variance,  a  is  the  corresponding 
power-law  exponent  of  the  fractional-frequency  noise  spectrum  Sy(f )  oc  fa.  In  the  context  here, 
its  valid  range  is  —4  <  a  <  2.  E  {TotHvar  (r,  T)}  relative  to  E  {Hvar  (r,  T)}  in  (6)  is  independent 
of  t  and  T,  dependent  on  noise  type,  and  biased  low,  giving  a  the  negative  sign  in  Table  1.  The  edf 
formula  (7)  is  a  convenient,  empirical  or  “fitted”  approximation  with  an  observed  error  below  10% 
of  numerically  computed  exact  values  derived  from  Monte-Carlo  simulation  method  using  the  bo  and 
b\  coefficients  of  Table  1  and  with  the  error  decreasing  with  averaging  factor  m  ~  t/tq  increasing. 
In  fact,  (7)  should  be  used  only  if  data-sampling  period  to  is  sufficiently  short  compared  to  the 
averaging  time  r  by  r/ro  >  16.  Otherwise,  there  are  not  enough  points  for  the  data-extension 
procedure  in  the  total  estimator  to  have  significant  advantage  over  the  plain  Hadamard  estimator. 
In  other  words,  the  ro-dependence  of  the  total  estimator  of  (4)  plays  a  significant  role,  whereas  the 
weaker  ro-dependence  of  the  maximum-overlap  estimator  of  plain  )  given  by  (2)  is  generally 

suppressed  as  in  (2).  It  is  well  known  that  maximum-overlap  statistical  estimators  will  increase 
edf,  hence  confidence,  and  the  degree  of  data  overlap  is  dependent  on  sampling  interval  to  relative 
to  t  [12, 13].  Real  data  should  be  sampled  as  fast  as  practical  for  a  given  averaging  time.  This 
is  especially  true  in  order  for  the  data  extension  of  each  subsequence  to  be  effective  in  the  total 
approach. 

Assuming  chi-square  distribution  properties  and  edf  computed  by  (7)  and  the  values  of  Table 
1,  confidence  intervals  will  be  conservative  since  the  distribution  is  actually  narrower  than  chi- 
square.  Although  not  quantitatively  investigated,  the  narrowing  of  the  distribution  is  proportional 
to  increasing  averaging  factor  m  —  t/to.  Fortunately  with  real  data  runs,  m  is,  of  course,  always 
largest  at  longest-term.  Depending  on  the  noise  type,  we  have  seen  narrowing  by  as  much  as  15% 
for  m  ~  100, 000. 

To  show  the  improvement  in  estimating  the  Hadamard  function,  Table  2  lists  the  exact  values 
of  edf  from  theory  for  computations  using  TotHvar  vs.  plain  Hvar  for  the  longest  averaging  factor  in 
which  t  =  T/3.  This  point  is  the  last  point  in  the  estimate,  and  the  improvement  in  confidence  using 


Table  1:  Coefficients  for  computing  (6)  and  (7),  normalized  bias  and  edf  of  TotHvar. 


Noise 

Abbrev. 

a 

a 

bo 

bi 

White  FM 

WHFM 

0 

-0.005 

0.559 

1.004 

Flicker  FM 

FLFM 

-1 

-0.149 

0.868 

1.140 

Random  Walk  FM 

RWFM 

-2 

-0.229 

0.938 

1.696 

Flicker  Walk  FM 

FWFM 

-3 

-0.283 

0.974 

2.554 

Random  Run  FM 

RRFM 

-4 

-0.321 

1.276 

3.149 

TotHvar  is  substantial,  particularly  for  the  general  case  of  WHFM  noise.  TotHvar  is  a  significantly 
improved  estimator  that  offsets  much  of  the  criticized  inefficiency  in  using  the  sample  Hadamard 
deviation  as  opposed  to  the  sample  Allan  deviation  in  the  presence  of  common  WHFM  noise  in 
frequency  standards. 
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Table  2:  Exact  gain  for  rmax  =  Tf  3. 


Noise 

edf  gain  of  TotHvar  (T/3,  T) 

WHFM 

3.447 

FLFM 

2.448 

RWFM 

2.044 

FWFM 

1.676 

RRFM 

1.313 

6  POWER  LAW  DETECTION 

It  is  important  to  be  able  to  determine  which  power-law  noise  type  is  present  for  a  given  r-value  in 
the  range  —  4  <  a  <  0  so  that  TotHvar’s  bias  can  be  removed  automatically.  Similarly,  before  the  edf 
can  be  determined  to  establish  confidence  intervals  and  set  error  bars  for  a  stability  measurement, 
it  is  necessary  to  identify  the  dominant  noise  process.  This  section  describes  a  noise-identification 
(noise-ID)  algorithm  that  has  been  found  effective  in  actual  practice,  and  that  works  for  a  single 
r-point  over  the  full  range  of  — 4  <  a  <  2.  It  is  based  on  the  Barnes  B1  function  [14],  which  is 
the  ratio  of  the  IV-sample  (standard)  variance  to  the  two-sample  (Allan)  variance,  supplemented 
by  applying  this  function  to  frequency  data,  and  the  R(n)  function  [15],  which  is  the  ratio  of  the 
modified  Allan  to  the  normal  Allan  variances. 

The  Bi  function  has  as  arguments  the  number  of  frequency  data  points  TV,  the  dead  time  ratio 
r  (which  is  set  to  1),  and  the  power-law  r-domain  exponent  /i.  The  Bl  dependence  on  p,  is  used 
to  determine  the  power-law  noise  type  for  — 2  <  p  <  2  (WHPM  and  FLPM  to  FWFM).  For  a  Bl 
that  is  consistent  with  a  p  =  — 2  result,  the  a  =  1  or  2  (FLPM  or  WHPM  noise)  ambiguity  can  be 
resolved  with  the  R(n)  ratio  using  the  modified  Allan  variance. 

For  the  Hadamard  variance,  the  noise  determination  must  be  extended  to  p  =  3  (or  a  =  —4, 
RRFM).  This  can  be  done  by  applying  the  Bl  ratio  to  frequency  (instead  of  the  usual  phase)  data 
and  adding  2  to  the  resulting  p.  This  procedure  is  called  “*B1”  herein.  Since  the  *B1  procedure 
simply  applies  the  Barnes  Bl  ratio  to  frequency  data  instead  of  phase  data,  its  use  is  as  before,  but 
now  its  range  is  effective  from  WHFM  to  RRFM  noise.  (This  is  analogous  to  simulation  of  RRFM 
data  by  treating  RWFM  phase  data  as  frequency  data.) 

The  overall  noise  identification  process  is  as  follows: 

•  calculate  the  standard  and  Allan  variances  for  the  applicable  r  averaging  factor, 

•  calculate  Bl,  B1(N,  r  =  1,  p)  = 

•  determine  the  expected  BI  ratios  for  a  =  -3  through  1  or  2, 

•  set  boundaries  between  them  and  find  the  best  power-law  noise  match, 

•  resolve  an  a  =  1  or  2  ambiguity  with  the  modified  Allan  variance  and  R(n),  or 

•  resolve  an  a  —  -3  or  -4  ambiguity  with  *B1. 


Table  3:  Formulae  for  131  (TV,  r  =  l,p).  Substituting  frequency  data  into  the  usual  phase-data 
measurement  of  Bl  ratio  will  shift  these  formulae  to  the  p  +  2  range,  thus  covering  RRFM. 
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Noise 

Bl  = 

FWFM 

2 

rnim+T) 

6 

RWFM 

1 

N/2 

FLFM 

0 

N  InN 
2(N-1 )  ln2 

WHFM 

-1 

1 

WH  or  FL  PM 

_o 

N2- 1 

Zi 

1.5(  Af)(JV— 1) 

For  a  data  run  of  length  N,  Table  3  gives  five  specific  formulae  for  B1  corresponding  to  p  = 
-2,  -1,  0,  1,  and  2.  Table  4  summarizes  the  power-law  detection  scheme  and  gives  the  boundaries 
for  demarcating  each  noise  type.  The  boundaries  between  the  Bl,  *B1,  and  R(n)  functions  are, 
in  general,  set  as  the  geometric  means  of  their  expected  values,  and  the  actual  measured  ratio  is 
tested  against  those  values  downward  from  the  largest  applicable  p.  For  example,  if,  during  the 
testing,  the  measured  Bl  ratio  is  greater  than  the  square  root  of  the  product  of  the  expected  Bl 
values  for  RWFM  and  FLFM  noise,  it  is  determined  to  be  the  former  ( a  =  —2,  RWFM). 

High  levels  of  frequency  drift  should  be  removed  to  best  identify  the  underlying  noise  process 
by  this  method.  Also,  the  R(n)  ratio  cannot,  of  course,  be  used  for  r  =  to  averaging  factor  (in 
which  case  it  is  1  for  all  noise  types).  Finally,  at  the  very  longest  averaging  factor  or  last  r-point, 
it  is  better  to  use  the  previous  or  r  —  ro  point  to  estimate  the  noise  type.  This  algorithm  has 
been  used  in  commercial  frequency-stability  software  [16]  for  the  past  decade  with  good  success. 
It  allows  bias  corrections  and  error  bars  to  be  calculated  automatically  during  an  analysis  for  all 
of  the  common  time-domain  stability  statistics  (including  the  new  Hadamard  total  variance  here) 
over  the  full  range  of  noise  types  and  for  essentially  all  r  averaging  times. 

7  THE  KALMAN  NOISE  MODEL  AND  THE  GPS  OPERA¬ 
TIONS  PROBLEM 

The  time  update  of  clock  states  in  the  Master  Control  Station  (MCS)  Kalman  prediction  algorithm 
is  based  on  an  average  of  the  the  most  recent  measurement  of  these  states  for  each  individual  clock, 
modeled  simply  by  random  noise  acting  on  phase  x(t),  frequency  y(t),  and  frequency  drift  z(t). 
With  this  model,  the  measured  power-law  a  exponents  of  the  frequency-fluctuation  noise  spectrum 
take  on  only  the  values  0,-2,  and  -4,  corresponding  to  WHFM,  RWFM,  and  RRFM,  or  p  =  -1, 
1,  and  3  in  the  r-domain.  Hence,  we  want  to  precisely  extract  the  level  of  these  noises  for  each 
clock  using  the  most  efficient  method  possible,  which  heretofore  has  been  the  sample  Allan  variance 
with  drift  removed  from  the  data  run,  and  more  recently  the  sample  Hadamard  variance,  because 
of  its  logical  link  to  the  model.  If  white  PM  (WHPM)  is  a  significant  noise  component,  and  for 
completeness,  the  a  =  2,  p  =  —  2  case  corresponding  to  WHPM  is  included  as  a  separate  error. 

The  parameters  used  by  the  MCS  within  GPS  system  operations  are  denoted  as  Kalman  filter 
q’s.  By  convention,  each  filter  parameter  qi,  i  =  0, 1,2,3  corresponds  respectively  to  r-domain 
power  law  exponents  p  =  —2,  —1, 1, 3.  For  the  Hadamard  variance,  the  relationship  is  [3] 

2/  \  2  2  2  2 
H CTy  (r)  =  VwHPM  +  aWHFM  +  a RWFM  +  ° RRFM 

=  ^Qor~2  +  qir-1  +  jlq2r  +  y^<?3  r3.  (8) 

For  the  Allan  variance,  the  relationship  is  [17] 

(jy2{T)  =  3qor~2  +  qir'1  +  | q2r  [+^q3r3]  ,  (9) 

where  the  inclusion  of  the  RRFM  noise  term  as  [+5U<?3T3]  is  a  point  of  contention  for  two  reasons. 
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Table  4:  Power-Law  Noise  Identification. 


Noise 

a 

A* 

ID  by 

Remarks 

RRFM 

-4 

3^ 

B1&*B1 

Use  *B1  to  resolve  a  =  -3  or  -4  ambiguity 

Decision  boundary:  (Bl(FWFM)  +  Bl(RWFM)}  /  2 

FWFM 

-3 

2 

B1&*B1 

Use  *B1  to  resolve  a  =  -3  or  -4  ambiguity 

Decision  boundary:  (Bl(FWFM)  +  Bl(RWFM)}  /  2 

1 

B1 

Decision  boundary:  sqrt  {Bl(RWFM)  x  Bl(FLFM)} 

-1 

B1 

Decision  boundary:  sqrt  (Bl(FLFM)  x  Bl(WHFM)} 

WHFM 

0 

mm 

B1 

Decision  boundary:  sqrt  (Bl(WHFM)  x  Bl(FLPM)} 

FLPM 

1 

-2 

Bl&R(n) 

Use  R(n)  to  resolve  a  =  1  or  2  ambiguity 

Decision  boundary:  sqrt  (Bl(FLPM)  x  Bl(WHPM)} 

WHPM 

2 

-2 

Bl&R(n) 

Noise  ID  Methods: 

B1  =  Barnes  B1(N,  r,  m)  bias  function  with  r  =  1  [14]. 

*B1  =  B1  applied  to  frequency  data  as  phase  data  with  q  =  p,  +  2.  R(n)  = 

ratio,  mod  Allan  variance/ Allan  variance.  [15]. 

First,  estimating  q%  by  (9)  using  real  data  is  unreliable  because  RRFM  is  inconsistent  by  the 
definition  of  the  Allan  variance.  Second,  ref.  [17],  from  which  the  term  derives,  does  not  compute 
the  Allan  variance;  instead,  it  computes  the  optimal  mean-square  prediction  error  variance  of 
y(to,to  +  t )  based  on  {x(t),t  <  to},  for  frequency  noise  spectra  with  a  =  0,  -2,  and  -4.  For  these 
reasons,  we  advise  omitting  the  RRFM  term  entirely  from  (9).  The  other  terms  of  (9)  happen  to 
be  correct  for  Allan  variance. 

The  GPS  Hadamard  variance  is  defined  to  be  equivalent  to  the  Allan  variance  for  WHFM, 
which  is  confirmed  in  comparing  (8)  and  (9);  however  the  variances  differ  by  a  factor  of  two  for 
RWFM,  therefore  they  cannot  be  used  interchangeably  under  normal  circumstances  and  involving 
drift-free  stochastic  processes. 

Tuning  the  Kalman  filter  depends  on  the  ability  to  “q”  each  individual  clock  according  to 
estimates  of  its  noise.  The  GPS  Block  HR  satellite  program  incorporates  Rb  atomic  oscillators 
that  are  characterized  by  a  mix  of  various  levels  and  types  of  random  noise  and  with  frequency 
drift  that  may  be  significantly  above  noise.  Tins  kind  of  oscillator  mix  is  difficult  to  manage  using 
Avar  and  (9),  which  must  be  used  based  on  drift-removed  frequency  residuals.  However,  reverting 
to  using  “frequency-drift  insensitive”  Hvar  and  using  (8),  the  confidence  becomes  a  factor  of  about 
g  less  near  the  last  and  crucial  long-term  rmax  =  T '/3  value  owing  to  the  plain  sample  Hadamard’s 
edf  of  one  less  as  compared  to  Allan’s  edf.  For  the  proper  perspective,  note  that  we  are  in  the 
one-week  averaging  r-region  with  a  last  real-time  data  run  of  about  one  month,  thus  edf  «  1- 
2;  so  estimating  filter  q' s  is  somewhat  subjective.  Figure  2  illustrates  a  summary  of  estimates  of 
frequency  stability  for  each  GPS  satellite  clock  as  published  in  reports  issued  by  the  Naval  Research 
Laboratory  [18]. 

Table  2  shows  that  the  new  TotHvar  (T/3,  T)  edf  is  multiplied  by  a  factor  of  1.3  to  3.4  over 
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plain  Hvar  (T/ 3,  T).  TotHvar  can  be  applied  directly  and  reliably,  while  retaining  the  efficiency  of 
the  sample  Allan  variance  without  the  difficulty  associated  with  real-time  drift  removal. 

The  work  of  this  paper  has  impact  on  two  GPS  operational  issues.  The  first  is  that  the  time 
needed  to  estimate  the  Hadamard  variance  is  substantially  reduced.  For  example,  to  obtain  a  r  = 
one-week  estimate  of  the  Hadamard  variance  with,  say,  the  last  40  days  of  measured  data,  the  Total 
approach  using  TotHvar  obtains  a  one-week  estimate  with  the  same  or  better  confidence  in  about 
26  to  34  days  of  measured  data.  The  second  issue  is  that  satellite  data  are  obtained  by  the  linked 
common- view  method  [19],  and  the  delay  in  receiving  the  monitor  station  tracking  data  is  currently 
at  2  to  3  days.  Thus,  it  is  important  to  extract  maximum  information  from  data  at  hand. 
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Figure  2:  Hadamard-deviation  frequency  stability  of  individual  GPS  satellite  clocks  vs.  USNO 
Master  Clock  for  the  period  1  January  to  1  July,  2000  [18]. 


8  EXAMPLE 

Figure  3  is  data  of  SV24,  a  Block  IIA  GPS  satellite.  Total  Hadamard  deviation,  plain  Hadamard 
deviation,  and  Allan  deviation  are  compared  with  increasing  data  spans  starting  at  7  days  and 
extending  to  28  days  and  shows  how  each  of  these  statistics  behaves  as  it  evolves.  As  is  generally 
the  case,  TotHdev  performs  better  at  estimating  the  longest-term  noise  level  than  plain  Hdev  for 
measured  data  spans  as  indicated  by  estimated  levels  from  later  (longer)  data  spans. 


9  CONCLUSION 

We  have  developed  a  significantly  improved  estimator  of  the  three-sample  Hadamard  frequency 
variance  based  on  the  so-called  “total”  approach  and  denoted  as  TotHvar,  for  use  in  GPS  operations 
and  analysis.  Practically  speaking,  we  have  reduced  the  long-term  estimation  uncertainty  in  terms 
of  edf  by  a  factor  of  1.3  to  3.4,  depending  on  the  noise  type,  and  we  have  presented  a  way  to 
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Figure  3:  Total  Hadamard  deviation,  plain  Hadamard  deviation,  and  Allan  deviation  for  SV24 
satellite  clock  data  as  the  data  run  increases  from  7  days  (front  plot)  to  28  days  (rear  plot).  The 
last  (rightmost)  values  of  TotHdev  for  shorter  data  runs  anticipates  the  underlying  noise  level  of 
longer  runs  compared  to  plain  Hdev  (arrowed  lines  are  projected  off  28-day  data  run).  The  Allan 
deviation’s  response  to  frequency  drift  masks  the  long-term  noise  level. 


automatically  remove  the  moderate  negative  bias  of  TotHvar  by  a  power-law  detection  algorithm. 
Having  confidence  greater  than  plain  Hvar  and  even  equal  to  or  greater  than  Avar,  TotHvar  is  a 
statistic  that  permits  tuning  of  the  MCS  Kalman  filter  with  more  accurately  chosen  clock-estimation 
parameters  (or  q’s)  that  are  linked  to  the  most  recent  measurements  of  frequency  stability  of 
each  clock.  The  increased  confidence  from  TotHvar  and  shorter  data  processing  delays  will  play 
significant  roles  in  adequately  managing  future  GPS  system  events. 
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Questions  and  Answers 


MASSIMO  TINTO  (JPL):  I  have  a  quick  question.  The  Hadamard  function,  in  a  sense,  can 
be  seen  as  the  Allan  variance  of  the  Allan  variance  because  you’re  taking  the  difference  of  the 
difference.  From  your  definition,  you  have  Y  (T)  minus  Y  (T)  plus  r,  minus  the  shift  of  1  by 
r.  So  have  you  thought  about  doing  the  Hadamard  of  the  Hadamard,  so  going  to  the  third 
order?  So  going  a  step  further  from  what  you  have  done  and  see  if  that  will  also  give  you 
some  extra— 

DAVID  HOWE:  I  haven’t.  The  problem  with  going  to  higher  differences  is  that  the  efficiency 
goes  down  for  the  kinds  of  noise  processes  which  are  characteristic  of  the  clocks,  because  you 
lose  degrees  of  freedom  as  you  have  to  use  longer  data  lengths  to  estimate  averaging  times. 

TINTO:  Oh,  I  see. 

STEVE  HUTSELL  (USNO  AMC):  Dave,  I  was  wondering  if  you  could  comment  on  the 
increase  of  the  effective  degrees  of  freedom.  How  much  of  it  was  due  to  the  overlapping 
technique  and  how  much  of  it  was  due  to  the  extension  in  the  total  technique? 

HOWE:  All  of  this  was  due  to  the  extension.  That  gain  was  due  to  the  extension  entirely 
because  the  numbers  were  taken  at  the  last  point.  There  is  not  overlapping  at  the  last  point. 

Now,  that  does  raise  the  point  that,  as  a  practical  matter,  total  and  overlapping  estimators  rely 
heavily  on  sampling  as  fast  as  practical.  You  understand  that  you  take  estimates  and  you  shift 
by  to>  your  sampling  interval.  Total  takes  advantage  of  that  more  fully  than  overlapping,  but 
that  has  been  true  for  quite  a  while. 
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